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Abstract 

Transient  thermal  analysis  plays  a  central  role  in  the  design  and  optimization  of  high  temperature  solid  oxide  fuel  cells  (SOFCs)  during 
startup/shutdown,  because  of  the  potential  for  damaging  thermal  gradients  to  develop  within  the  SOFC  components.  The  optimal  design  of  a 
heating/cooling  process  is  one  that  minimizes  the  total  time  required  to  reach  a  prescribed  final  operating  temperature,  while  not  exceeding 
given  thresholds  of  maximum  allowable  temperature  gradients.  To  this  end,  we  consider  the  SOFC  unit  cell,  which  is  heated  by  hot  air  supplied 
into  the  oxidizer  channel  at  a  specified,  time-dependent  inlet  temperature.  Beginning  with  a  general  thermal  model  of  the  cell,  we  develop  and 
evaluate  limiting  cases  that  allow  closed-form  analytical  solutions  of  the  time-varying  temperature  fields,  from  which  heating  time  and  maximum 
temperature  gradient  are  calculated.  The  results  are  generalized  by  presentation  in  terms  of  a  modified  effective  Peclet  number  and  dimensionless 
inlet  temperature  function.  Finally,  the  accuracy  of  these  predictions  is  evaluated  by  comparison  to  results  of  3-D  CFD  modeling  in  Fluent,  and 
design  maps  for  optimizing  the  transient  heating  process  are  presented.  Results  indicate  that  the  reduced-order  models’  simplicity,  computational 
savings,  and  ability  to  capture  the  essential  physics  of  the  transient  process  justify  their  use  in  design  calculations  over  more  complex,  highly 
detailed,  numerical/CFD  schemes. 

©  2005  Elsevier  B.V.  All  rights  reserved. 

Keywords:  Solid  oxide  fuel  cell;  Thermal  modeling;  Transient  analysis 


1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  are  currently  being  developed 
for  mobile  and  stationary  power  plant  applications,  and  much 
attention  is  being  paid  to  the  design  and  optimization  of  their 
steady  state  performance  in  an  effort  to  make  SOFC  technol¬ 
ogy  commercially  viable  in  the  near  future.  However,  it  is  being 
increasingly  realized  [1]  that  because  of  substantial  thermome¬ 
chanical  stresses  developed  in  the  cell  components  and  stack, 
the  transient  process  of  heating  an  SOFC  from  room  temper¬ 
ature  to  operating  temperature  (600-800  °C)  during  startup,  or 
cooling  down  to  ambient  during  shutdown  must  be  given  spe¬ 
cial  attention  as  well.  The  ultra  thin  electrolyte  and  electrode 
layers  (PEN  structure)  are  prone  to  delamination,  cracking,  or 
other  catastrophic  failure  if  subjected  to  excessive  thermal  shock, 
temperature  gradients,  and  thermal  cycling  during  startup  or 
shutdown  processes. 
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These  dangers  can  usually  be  avoided  altogether  by  proceed¬ 
ing  through  the  transient  in  a  very  slow,  controlled  fashion  (in  a 
recent  ASME  conference  presentation  Hawkes  et  al.  [2]  reported 
taking  2  days  to  bring  a  stack  up  to  operating  temperature).  As 
SOFC  technology  matures,  however,  it  is  likely  that  the  con¬ 
sumer  will  demand  that  the  fuel  cell  reach  operating  conditions 
as  quickly  as  possible,  for  example,  in  the  cold  start  of  an  auto¬ 
mobile  [3].  Thus,  the  optimal  design  of  a  startup  process  will 
minimize  the  total  time  required  for  heating,  subject  to  the  con¬ 
straint  of  some  maximum  allowable  temperature  gradient  (to 
avoid  thermomechanical  fracture)  as  well  as  time  rate  of  change 
of  temperature  (to  avoid  thermal  shock  and  creep  stresses).  It 
will  also  be  necessary  to  quantify  the  expected  number  of  ther¬ 
mal  cycles  the  cell  can  withstand  during  its  lifetime,  subject  to 
the  real  possibility  that  faster  heating  times  will  lower  the  life 
expectancy  of  the  cell. 

In  the  literature,  some  preliminary  efforts  to  address  these 
concerns  have  been  reported.  Although  numerous  papers  have 
presented  the  development  of  numerical/CFD  models  to  simu¬ 
late  SOFC  behavior  at  steady  state,  few  are  capable  of  simu¬ 
lating  transient  operation  [1,3-5],  and  only  one  [1]  has  begun 
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Nomenclature 

A  cross-sectional  area  of  a  cell  component  (m2) 

cp  specific  heat  (J  kg- 1  K- 1 ) 
h  convective  heat  transfer  coefficient  (Wm-2  K-1) 

k  thermal  conductivity  (Wm-1  K-1) 

K  rate  of  inlet  temperature  rise  (°C  s-1) 

Kq ff  dimensionless  (effective)  rate  of  inlet  temperature 

rise 

L  length  of  channel  (m) 

Pe  effective  Peclet  number  of  the  flow,  uQffL/aQff 

7f  final  (operating)  temperature  (°C) 

To  initial  temperature  (°C) 

u  velocity  of  hot  air  stream  (ms-1) 

ue  ff  effective  velocity  of  hot  air  stream  (ms-1) 

Greek  letters 

a  thermal  diffusivity  (m2  s-1) 

c^eff  effective  thermal  diffusivity  (m2  s-1) 

y0eff  parameter  relating  conduction  to  advection  (m) 

r\  kinematic  viscosity  (m2  s-1) 

p  density  (kg  m- 3 ) 

rc  advective  time  scale  (s) 

Th  total  heating  time  (s) 

dimensionless  total  heating  time 
Ti  inlet  temperature  function  time  scale  (s) 


to  quantify  thermally  induced,  transient  thermal  stresses  during 
startup/shutdown.  These  models  are  based  on  finite  element  or 
volume  approaches  with  the  intent  to  give  highly  detailed  infor¬ 
mation  about  current  density,  species  distribution,  flow  and  tem¬ 
perature  fields,  and  propagation  of  thermal  waves  in  the  solid. 
Thus,  the  information  required  for  transient  design  would  be 
available,  but  the  complexity  of  the  models  and  the  copiousness 
of  the  results  may  obscure  the  underlying  physics  of  the  process 
and  prohibit  the  development  of  simple  design  rules  for  opti¬ 
mal  transient  heating/cooling  processes  (assuming  such  rules  do 
exist).  In  addition,  the  rigorous  CFD-based  approaches  demand 
a  great  deal  of  computational  power,  especially  when  simulat¬ 
ing  a  long  transient  process  with  a  large  number  of  parameters 
to  study.  Adding  to  the  computational  expense  is  the  very  fine 
discretization  required  for  numerical  modeling  because  of  the 
very  high  aspect  ratio  of  components  in  the  cell,  for  example, 
the  15  [xm  thick  by  10  cm  long  electrolyte  layer. 

To  overcome  the  above-mentioned  conceptual  and  compu¬ 
tational  problems,  we  aim  at  developing  reduced-order  models 
and  an  analytical  approach  leading  to  the  closed-form  solution 
of  the  problem.  Our  ultimate  goal  is  to  develop  simple,  efficient 
design  rules  that  clearly  show  the  effects  of  the  various  system 
parameters  on  (1)  the  total  time  required  for  startup  and  (2)  the 
maximum  temperature  gradients  developed  during  the  process. 
The  interest  is  not  in  developing  a  model  that  can  give  a  highly 
detailed  prediction  of  the  temperature  field  at  any  given  moment 
in  time,  but  rather,  a  model  that  accurately  and  efficiently  pre¬ 
dicts  the  global  quantities  just  mentioned.  This  paper  is  the  first 


to  attempt  to  develop  an  efficient  design  strategy  for  optimizing 
the  transient  process  and  will  focus  on  the  initial  startup  of  a 
cell  from  ambient  temperatures  with  a  simplified  approach  that 
generally  (with  few  modifications)  applies  to  shutdown  as  well. 

2.  Model  formulation  and  approach 

Under  consideration  is  the  planar  type  SOFC,  which  is  a  stack 
of  repeating  unit  cells  with  dimensions  shown  in  Fig.  1.  The 
interconnects  (current  collectors)  are  made  of  stainless  steel,  the 
solid  electrolyte  is  yttria- stabilized  zirconia  (YSZ),  the  porous 
anode  is  Ni-doped  YSZ,  and  cathode  is  Sr-doped  LaMnC>3 
(LSM).  The  cell  is  heated  by  flowing  hot  air  into  the  oxidizer 
channel  while  controlling  the  inlet  air  temperature  as  a  function 
of  time.  As  a  first  approximation,  conditions  in  the  fuel  channel 
are  assumed  quiescent  (negligible  flow),  with  the  composition 
of  the  gas  phase  similar  to  what  would  be  found  in  a  typical  fuel 
stream.  At  startup,  no  electrochemical  reactions  take  place  until  a 
prescribed  temperature  is  reached,  at  which  point  electrochemi¬ 
cal  “light-off”  occurs  [6].  The  electrochemical  process  generates 
heat,  and  consumes  and  produces  chemical  species  in  the  flow 
streams,  coupling  the  thermal-fluid  model  to  the  electrochemical 
model.  For  simplicity,  it  is  assumed  that  a  desired  operating  tem¬ 
perature  is  reached  without  triggering  electrochemical  reactions 
(i.e.  no  heat  generation),  so  that  attention  is  limited  to  convective 
heating  of  the  cell. 

Another  simplifying  assumption  is  that  of  a  1-D  tempera¬ 
ture  field  in  each  component  of  the  cell,  with  variation  only 
along  the  flow  direction,  as  suggested  by  the  high  aspect  ratio 
of  the  channel  (>30:1).  The  small  physical  dimensions  of  the 
components  also  suggest  that  they  may  be  in  local  thermal 
equilibrium  (normal  to  the  flow  direction)  and  this  possibil¬ 
ity  is  investigated  through  the  development  and  analysis  of 
three  reduced-order  models  of  varying  complexity.  The  first  and 
most  general  model  assumes  the  solid  temperature  is  different 
from  the  gas  temperature,  and  the  resulting  two-equation,  cou¬ 
pled,  transient  model  is  solved  numerically.  The  second  model 
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Fig.  1.  Geometry  of  the  unit  cell  of  a  planar  type  SOFC  stack. 
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Fig.  2.  Schematic  of  the  unit  cell  as  a  channel  with  composite,  insulated  walls. 
The  temperature  profile  is  assumed  1-D  in  each  layer  with  variation  in  the  flow, 
z-direction  only. 


where  h  is  the  convective  heat  transfer  coefficient,  Tq  and  T\c 
are  temperatures  of  the  cathode  and  interconnect,  respectively, 
Pg-c  and  Pg_ic  are  the  contact  perimeters  between  the  air  (gas) 
channel  and  the  cathode  layer,  and  the  air  (gas)  channel  and  the 
interconnect  layer,  respectively,  A  is  the  cross  sectional  area  of 
the  air  channel,  and  k ,  p,  cp  are  the  thermal  conductivity,  density, 
and  specific  heat  of  the  air,  respectively.  Similar  equations  are 
written  for  the  other  layers,  for  example,  the  interconnect  layer 


(pep  A)j^ 


dTlc 

dt 


02  j1 

(*A)IC-^+APg-ic(rg 

Pic-c 


Tic) 


Ric-c 


( Tic  ~  TC) 


(3) 


which  is  coupled  to  the  air  equation  through  the  convection  term 
(the  second  term  on  the  right-hand-side),  and  to  the  cathode  layer 
(see  Figs.  1  and  2)  through  the  thermal  resistance 


assumes  the  gas  and  solid  are  in  local  thermal  equilibrium  result¬ 
ing  in  a  single  transient  convective-conductive  equation  for 
which  an  analytical  solution  has  been  obtained  [7].  The  third 
model  goes  one  step  further  in  simplifying  the  problem  and 
assumes  that  conduction  along  the  flow  direction  is  negligible 
compared  to  advection  of  thermal  energy  down  the  length  of 
the  cell.  This  latter  model  is  thus  a  purely  convective  heating 
model,  also  yielding  an  analytical  solution.  The  range  of  validity 
of  these  models  is  then  established  by  comparing  their  predic¬ 
tions  of  heating  time  and  maximum  temperature  gradients  to 
the  results  of  detailed,  3-D  CFD  simulations  of  SOFC  unit  cell. 
It  is  assumed  that  the  geometry  of  the  cell,  material  selection, 
and  initial  and  final  temperatures  are  prescribed  design  param¬ 
eters,  leaving  the  inlet  air  temperature  and  air  velocity  (flow 
rate)  as  parameters  that  can  be  used  to  optimize  the  heating 
process. 

2.7.  Two-equation,  coupled,  solid-gas  model 

The  most  general,  reduced-order,  transient  heating  model  we 
developed  is  derived  by  applying  conservation  of  energy  to  each 
layer  (component)  of  the  cell  shown  in  Fig.  2.  This  yields  a  set 
of  seven  partial  differential  equations  (PDEs)  that  are  coupled 
through  the  temperature  difference  between  adjoining  layers. 
The  equation  for  the  air  channel,  assuming  constant  velocity 
plug  flow  and  neglecting  thermal  radiation,  is 


(pcvA) 


377  37) 

—  +  u 
dt 


dz 


dzT0 


=  (kA) -  hPg-ciTg  -  7c)  -  hPg-iciTg  -  Tic )  (1) 


subject  to  the  boundary  and  initial  conditions 

dT 

BC’s  :  Tg(0,  t)  =  f(t);  — *(L,f)  =  0, 

oz 

Tg(z,  0 )  =  T 


IC 


(2) 


he  tc 
Ric-c  =  + 


2kic  Tkc 


(4) 


where,  T\ c,  tc  are  the  thickness  of  the  interconnect  and  cathode 
layers,  respectively.  In  like  manner,  the  remaining  equations  for 
the  cathode  (C),  electrolyte  (E),  anode  (A),  fuel  channel  (f),  and 
lower  interconnect  (IC2),  are 


(pCpA)c^  =  (M)c5  +  hPg-c(Ts  -  TC) 


dt 


+ 


dz2 

PlCi-C 


dTE 

(pcvA\—  =  (M)e 


^IC-C 

d2TE 


(Tic  -  Tc)  -  T^(Tc  -  TE), 


Rc-e 


dt 


+ 


dz2 

Pc-E 


Rc-e 


(TC  -  Te)  +  ^^(Te  -  7a), 
TvE— A 


37a  327a  Pe-a 

( pcpA)A =  (kA) +  — — —  (Te  -  7a) 
dt  dzz  Re-A 

-  hPf-A(T{  -  7X)  -  ^-^2- (7a  -  7ic2), 

TvA— IC2 

37fuei  32  7fuei 

(^PA)fuel^  =  (kA)fud^P  +  hPf-A(Tf  ~  7a) 


dt 


dz2 

+  h  Pf-ic2  ( Tf  —  Tic2  ) , 


dTic?  /T  ,,  92Tic9  ,  Pa-ic2 


(pcv  A\Ci  =  (kA)lc 


+ 


2  dz2  '  Ra-ic2 

—  hPf-ic2(Tf  -  Tic2) 


(Ta  ~  Tic,) 


(5) 


The  boundary  and  initial  conditions  for  Eqs.  (3)  and  (5)  are  all 
of  the  same  form 


BC  s 


37-  37) 

\0,t)=  -±(L,t)  =  0: 


dz 


dz 


IC  :  7i(z,  0)  =  To 

(6) 
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Table  1 


Material  properties  of  cell  components 


Component 

P  (kgm  3) 

k  (W  m-1  K-1) 

cpakg-'K-1) 

Anode 

3030 

5.84 

595 

(Ni-doped  YSZ) 

Cathode  (LSM) 

3310 

1.86 

573 

Electrolyte  (YSZ) 

5160 

2.16 

606 

Current  collector 

8030 

20.0 

502 

(SS) 

Air  channel 

0.58 

0.047 

1051 

Fuel  channel 

0.2 

0.2 

5000 

This  set  of  coupled  PDEs  can  be  solved  numerically  in  its  cur¬ 
rent  form,  but  a  reasonable  simplifying  assumption  at  this  point, 
is  that  the  temperatures  in  the  solid  are  locally  uniform  at  each 
cross-section  in  the  direction  normal  to  the  flow.  This  is  verified 
by  calculating  the  Biot  number  Bi,  for  this  configuration  as  the 
ratio  of  thermal  resistance  across  the  entire  composite  channel 
wall,  to  convective  thermal  resistance  between  the  channel  walls 
and  air  stream 


Bi  = 


^conduction 


R 


convection 


where  K  is  the  rate  of  temperature  rise  at  the  inlet  in  units  of 
°C  s_1 .  It  should  be  noted  that  although  the  results  for  the  linear 
temperature  increase  are  discussed  in  this  paper,  all  mathemat¬ 
ical  developments  are  general  and  could  be  readily  extended  in 
a  straight-forward  fashion  to  investigate  any  functional  depen¬ 
dence  of  the  inlet  air  temperature. 

Eqs.  (8)  and  (10)  are  solved  numerically  using  a  Crank- 
Nicholson,  finite  difference  scheme  [8].  The  gas  equation  for 
the  next  future  time  step  is  solved  by  guessing  a  temperature 
distribution  in  the  solid,  and  using  the  tri-diagonal  matrix  algo¬ 
rithm  [8]  to  invert  the  coefficient  matrix.  Using  this  temperature 
distribution,  the  solid  equation  is  solved  in  a  similar  manner,  and 
the  process  repeated  until  the  temperatures  converge  before  pro¬ 
ceeding  to  the  next  time  step.  The  calculations  proceed  forward 
in  time  until  steady- state  is  reached.  Temperature  gradients  (spa¬ 
tial  and  temporal)  are  calculated  numerically  from  the  resulting 
temperature  field  history,  allowing  the  maximum  temperature 
gradients  developed  through  the  entire  heating  process  to  be 
identified.  The  only  modification  necessary  for  simulating  cool¬ 
ing  (rather  than  heating)  processes  is  to  use  negative  values  for 
K  in  the  boundary  condition,  Eq.  (1 1). 

2.2.  Convective-conductive  model 


where  h  is  the  convection  heat  transfer  coefficient,  t  is  the  thick¬ 
ness  of  layer  i  and  k  is  thermal  conductivity.  Using  values  found 
in  Table  1  for  ku  dimensions  from  Fig.  1  for  tu  and  estimating 
60  Wm-2  K-1  for  fully-developed  internal,  laminar  air  flow 
in  a  3  mm  pipe,  yields  Bi  ~  0.02,  validating  the  assumption  that 
the  solid  can  be  treated  as  a  locally  isothermal  lumped  capaci¬ 
tance  element.  This  understanding  allows  us  to  combine  Eqs.  (3) 
and  (5)  to  formulate  a  single  energy  equation  for  the  effective 
temperature  of  the  solid  matrix 

y>cpA)^  =  Y,{kA)‘^f + hpgs(F  - T s)  (8) 

i  i 


with  boundary  and  initial  conditions 


BC’s 


dTs 

~7r~(L,  t)  =  0; 

oz 


IC  :  Ts(z,  0)  =  T0 

(9) 


The  equation  for  the  gas  is  similar  to  Eq.  (1) 


(pcpA) 


dTj 

dt 


u 


97g 

dz 


d2T0 


=  (kA)g^  +  hPg-s(Ts 


To  further  simplify  Eqs.  (8)  and  (10),  and  to  enable  closed- 
form  analytical  solution  of  the  problem,  an  assumption  of  local 
thermal  equilibrium  between  the  gas  and  solid  is  employed.  This 
bold  simplification  is  applied  to  the  problem  at  hand,  not  to  prove 
or  disprove  that  the  solid  and  gas  are  at  the  same  temperature,  but 
to  determine  whether  the  approach  will  yield  accurate  predic¬ 
tions  of  heating  time  and  maximum  temperature  gradients.  If  the 
model  can  accomplish  this  task,  and  do  it  analytically,  then  it  is 
a  very  powerful  design  tool  for  optimizing  the  transient  heating 
and  cooling  process. 

To  implement  the  local  thermal  equilibrium  assumption  cor¬ 
rectly  (i.e.,  without  even  locally  violating  energy  conservation), 
the  following  procedure  is  employed:  First,  both  conservation 
Eqs.  (8)  and  (10)  are  added  together,  which  cancels  the  cou¬ 
pling  convective  solid-gas  heat  exchange  term.  Then,  by  defi¬ 
nition  of  local  thermal  equilibrium,  the  temperatures  of  the  gas 
and  solid  are  made  the  same,  leading  to  the  following  single 
convective-conductive  governing  equation  describing  transient 
heating  dynamics  of  the  unit  cell 


dT 

~dt 


+  Weff 


dT 

dz 


d2T 

“eff9? 


(12) 


where,  weff  is  the  effective  velocity  and  aeff  is  the  effective  thermal 
diffusivity  defined  as  follows: 


with  the  same  boundary  and  initial  conditions  given  in  Eq.  (2). 
We  further  consider  the  case  when  the  inlet  air  temperature  is 
defined  by  a  linear  increase  in  time  until  a  desired  final  (Tf) 
steady- state  temperature  is  reached 


To  +  Kt 


t  < 

t  > 


(Tf  -  T0) 


K 

( Tf  ~  To) 


(11) 


(pcpA)g  ...  „  _  «• 

^  cff  —  r T  W,  O'cff  —  ^ — ,  (13) 

yypcpA),.  2J^cpA); 

i  i 

Note  that  the  summations  are  now  indexed  to  include  every  com¬ 
ponent  of  the  cell  including  flow  channels.  The  effective  velocity 
is  the  physical  air  velocity  scaled  by  the  ratio  of  heat  capacity 
of  the  air  channel  to  energy  stored  in  the  channel  walls.  The 
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effective  thermal  diffusivity  is  the  cross-sectional  area- weighted 
diffusivity  of  the  composite  channel  wall.  Both  of  these  terms 
are  dominated  by  the  relatively  massive  interconnects  with  their 
high  conductivity,  density,  and  cross-sectional  area. 

The  same  procedure  used  to  derive  Eq.  (12)  can  be  applied  to 
the  boundary  and  initial  conditions  (they  must  be  first  expressed 
in  a  conservative  basis  in  terms  of  the  total  energy  rate/flux  [7]), 
yielding 

3T 

EC’s  :  7X0,  t)  -  Aar— (0,  t)  =  /(/); 

az 

3T 

— (L,t)  =  0,  IC:  T(z,  0)  =  T0  (14) 

oz 

where  the  coefficient  /3eff  in  Eq.  (14)  combines  thermal  energy 
conduction  in  the  solid  to  thermal  energy  transport  in  the  air 
channel 


(^A)iCi  +  (^)c  +  (kA)E  +  (*A)a  +  (kA)f  +  (kA)  ic2 


(pcvA)gu 


(15) 


Notice  that  if  heat  diffusion  in  the  air  is  neglected  -  a  very 
reasonable  approximation  as  compared  to  advection  -  then, 
fieff  =  OiQ  ff/ UQff  • 

When  temperature,  length,  and  time  are  normalized  using  the 
following  characteristic  scales: 


T-Tp 

Tf-To' 


(L  /  Uq  ff) 


(16) 


the  dimensionless  form  of  Eq.  (12)  can  be  expressed  as 


ar*  ar*  1  d2r* 

dt*  +  dz*  Pe  dz*2 


(17) 


where  Pe  is  defined  as  the  effective  Peclet  number, 
Pe  =  weff  L/aeff,  that  the  rati°  of  advection  of  thermal  energy 
down  the  length  of  the  channel  to  diffusion  of  thermal  energy  in 
the  solid  layers  of  the  SOFC  unit  cell.  Likewise,  the  boundary 
and  initial  conditions  are 


1  dT* 

BC’s  :  T*( 0,  t*)  -  —  —  (0,  f)  =  F(t*); 


dT 


* 


(U*)  =  0, 


3  z 


* 


Pe  dz* 

IC:  T*(z*,0)  =  0 


(18) 


where  F(t  )  is  the  non-dimensional  form  of  the  inlet  temperature 
function  given  by  Eq.  (11) 


F(f)  = 


KL 


(7f  -  7o)weff 


(19) 


Eq.  (19)  reveals  the  second  dimensionless  parameter  that  deter¬ 
mines  the  thermal  response  of  the  cell,  the  effective  rate  of  inlet 
temperature  rise 


E/ Uq  ff 

(7f  -  T0)/K 


(20) 


Physically,  this  parameter  represents  the  ratio  of  the  advective 
time  scale,  rc  =  L/weff,  to  the  time  scale  associated  with  the  tran¬ 
sient  forcing  of  the  inlet  temperature,  x\  -  (7f  —  Tq)/K  that  is, 


the  time  required  for  the  inlet  air  temperature  to  reach  the  final 
desired  temperature.  A  significant  insight  into  the  transient  heat¬ 
ing  process  can  be  obtained  by  analyzing  this  parameter.  In 
advection  dominated  flows  ( Pefy>  1,  note  the  modified  defini¬ 
tion  of  effective  Peclet  number  here,  combining  properties  of 
the  air  and  solid),  a  value  of  KQ ff  greater  than  unity  implies  that 
the  rate  of  thermal  energy  input  at  the  inlet  of  the  cell  exceeds 
the  capability  of  the  cell  to  store  and  distribute  this  heat  input. 
Thus,  temperature  at  the  inlet  builds  up  too  quickly,  leading  to 
excessive  thermal  gradients  without  significant  improvement  in 
total  heating  time. 

The  usefulness  of  Eq.  (17)  lies  in  the  fact  that  an  analyt¬ 
ical  series  solution  has  been  found  [7]  (see  Appendix  A  for 
details),  and  that  the  temperature  distribution  dependence  has 
been  reduced  to  only  two  dimensionless  parameters:  effective 
Peclet  number,  Pe ,  and  the  effective  rate  of  inlet  temperature 
rise,  Kq ff.  Also,  this  solution  provides  a  convenient  way  to  val¬ 
idate  the  numerical  scheme  used  to  solve  Eqs.  (8)  and  (10)  in 
Section  2. 1 .  As  the  heat  transfer  coefficient  coupling  the  gas  and 
solid  temperatures  is  made  arbitrarily  large,  the  temperatures  of 
the  gas  and  solid  should  converge  onto  each  other  and  the  result¬ 
ing  single  temperature  solution  must  match  the  analytical  result 
given  by  the  closed-form  solution  of  Eq.  (17). 


2.3.  Purely  convective  model 

A  very  interesting  limiting  case  of  Eq.  (17)  is  that  of  purely 
convective  flow  (Pe  — >  oo),  in  which  case  the  right-hand  side  of 
Eq.  (17)  vanishes.  The  resulting  equation 


ar*  (  ar* 

dt*  +  3z* 


and  in  particular,  its  dimensional  form 


dr 

~dt 


+  Weff 


dr 

dz 


(21) 


(22) 


looks  very  much  like  the  governing  equation  for  a  cell  with 
thermally  thin  cell  components  (i.e.  no  thermal  energy  storage 
in  cell  components,  and  local  thermal  equilibrium  between  the 
solid  and  gas  [7]).  However,  because  the  effective ,  not  physical 
velocity  is  used  here,  thermal  energy  storage  in  the  cell  is  prop¬ 
erly  accounted  for,  and  only  thermal  energy  conduction  along 
the  axial  direction  is  neglected.  The  relevant  boundary  and  initial 
conditions  are,  respectively 


r*(o,f*)  =  F(t*y,  r*(z*,o)  =  o 


(23) 


The  method  of  characteristics  yields  the  analytical  solution  of 
the  time  dependent,  1-D  temperature  distribution 


T*(z*,t*) 


0,  z*  >  t* 

KQff(t*  -  z*),  z*  <  t* 


(24) 
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This  very  simple  equation  for  the  unsteady  temperature  distribu¬ 
tion  provides  simple  algebraic  relationships  for  dimensionless 
heating  time 


1 


K,  ff 


+  1 


and  dimensionless  temperature  gradients 


(25) 


max 


dr* 


dt 


* 


=  r  eff 


max 


(26) 


The  dimensional  form  of  these  equations  yields  some  physical 
insight.  The  total  heating  time  is  the  sum  of  the  two  time  scales, 
rc  and  Ti  previously  discussed  in  Section  2.2  above 


rh 


Tf-Tp  |  L 

K  Uq  ff 


(27) 


The  fastest  that  the  cell  can  be  heated,  assuming  it  is  a  per¬ 
fect  heat  exchanger,  is  given  by  the  time  required  to  bring  the 
incoming  air  to  the  desired  operating  temperature,  plus  the  time 
required  for  thermal  energy  to  travel  from  the  inlet  to  the  exit. 
So,  Eq.  (27)  is  the  theoretical  lower  bound  (minimum  possible 
value)  on  heating  time,  and  unambiguously  suggests  that  heating 
time  is  inversely  proportional  to  both  K  and  u.  The  maximum 
dimensional  (physical)  temperature  gradients  become 


dT 

dz 


max 


K 

U  eff 


dT 

~dt 


=  K 

max 


(28) 


indicating  that  spatial  temperature  gradient  is  inversely  propor¬ 
tional  to  u ,  but  both  temperature  gradients  (spatial  and  temporal) 
increase  linearly  with  K.  Eq.  (28)  set  upper  bounds  on  the  max¬ 
imum  temperature  gradients  that  can  exist  in  the  solid  during 
the  transient  heating,  owing  to  the  fact  that  thermal  conduc¬ 
tion,  which  tends  to  diminish  temperature  gradients,  has  been 
neglected  in  the  formulation. 

Eqs.  (27)  and  (28)  constitute  very  simple  design  rules  that 
establish  relationships  between  design  parameters  and  clearly 
explain  the  competing  effects  that  must  be  balanced  to  optimize 
the  heating  process 


•  Increasing  flow  velocity,  u,  tends  to  decrease  both  the  heating 
time  ( positive  effect )  and  spatial  temperature  gradient  ( posi¬ 
tive  effect ),  but  has  no  effect  on  the  temporal  gradient. 

•  Increasing  the  rate  of  inlet  temperature  rise,  K ,  tends  to 
decrease  the  required  heating  time  ( positive  effect ),  but 
increases  both  the  temporal  and  spatial  temperature  gradients 
( negative  effect ). 


Thus,  the  purely  convective  model  has  yielded  useful  infor¬ 
mation  about  limiting  cases,  which  in  hindsight  seems  to  be 
almost  intuitively  obvious,  and  provides  the  framework  for  opti¬ 
mization  of  the  heating  process,  even  though  the  validity  of 
the  results  still  needs  to  be  established.  This  is  done  in  the 
next  section.  In  order  to  determine  whether  or  not  the  sim¬ 
ple,  reduced-order  transient  heating  models  have  any  significant 
bearing  on  reality,  we  compared  their  predictions  of  heating  time 


and  temperature  gradients  to  the  results  of  transient  3-D  CFD 
simulations  of  SOFC  unit  cell  heating. 

2 A.  CFD  model 

To  critically  compare  the  predictive  power  of  the  reduced- 
order  transient  thermal  models,  a  fully  three-dimensional 
thermal-fluid  analysis  of  the  unit  cell  in  Fig.  1  was  performed 
via  the  commercially  available  CFD  software  Fluent.  In  simu¬ 
lations,  thermophysical  properties  of  the  solid  materials  were 
assumed  constant  (estimated  at  the  average  temperature)  and 
are  listed  in  Table  1.  The  steady-state  flow  field  in  the  oxidizer 
channel  was  found  under  the  assumption  of  laminar  flow  using 
constant  thermophysical  fluid  properties  evaluated  at  the  average 
temperature  of  300  °C.  For  simplicity,  this  flow  field  was  then 
used  for  the  unsteady  calculations  of  the  temperature  field  in 
the  cell,  as  the  inlet  temperature  was  linearly  increased  from  25 
to  625  °C.  Once  the  inlet  air  temperature  reached  625  °C,  it  was 
held  constant  until  the  normalized  temperature  of  the  solid  at  the 
end  of  the  flow  channel  was  within  5%  of  the  normalized  steady 
state  value  (i.e.  595  °C),  at  which  point  the  simulations  were  ter¬ 
minated.  The  temperature  history  of  the  solid  (electrolyte)  was 
recorded  throughout,  and  post-processed  to  yield  maximum  spa¬ 
tial  and  temporal  temperature  gradients.  Because  of  symmetry 
along  the  axial  midplane,  only  half  of  the  unit  cell  was  mod¬ 
eled  using  44,000  discrete  volumes.  The  mesh  and  time-step 
were  properly  refined  to  ensure  that  the  results  were  mesh  and 
time- step  independent.  The  time- step  for  each  simulation  was 
adjusted  based  on  the  rate  of  inlet  temperature  rise,  K,  and  typ¬ 
ically  corresponded  to  a  2  °C  (per  time  step)  temperature  rise  at 
the  inlet. 

Cases  were  run  for  lvalues  of  0.01,  0.02,  0.05,  0.1,  0.2,  0.5, 
1.0,  and  2.0  °C  s_1 .  The  inlet  velocities  used  were,  1,  5,  10,  and 
20  m  s-1 ,  corresponding  to  effective  Peclet  numbers  of  0.8,  4.2, 
8.4,  and  16.7.  Results  of  these  simulations  are  compared  with 
results  of  the  reduced-order  models  in  the  next  section. 


3.  Model  results  and  analysis 


First,  the  results  of  CFD  simulations  were  compared  to  pre¬ 
dictions  obtained  with  the  most  general  reduced-order  model, 
which  is  the  two-equation,  coupled,  solid-gas  model  described 
in  Section  2.1.  Before  the  two-equation  model  can  be  numeri¬ 
cally  solved,  the  heat  transfer  coefficient,  h ,  must  be  specified 
at  each  point  along  the  air  channel.  Because  almost  one-third  of 
the  channel  is  in  the  thermally-developing  entrance  region,  the 
correlation  based  on  the  Graetz  solution  for  thermally  develop¬ 
ing,  hydrodynamically  fully-developed  laminar  flow  in  a  duct 
with  constant  wall  temperature  [9]  was  used  to  approximately 
account  for  variation  in  heat  transfer  coefficient  along  the  chan¬ 
nel  length 


Nux  ~  3.66  + 


0.0018 

v*1/3(0.04  +  v*2/3)2 


*  =  x/Dh 
Re  Pr 


(29) 


where  the  Nusselt  number  at  a  given  position,  x,  is  Nux  =  hxD\Jk , 
hydraulic  diameter  is  D\l  =  4Ac/Pw ,  the  Reynolds  number  is 
Re  =  uD\Jr\  and  Prandtl  number  is  Pr  —  rj/a. 
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Analytical 

Numerical 


Flow  Direction,  m 

Fig.  3.  Validation  of  numerical  solution  of  two-equation,  coupled  solid-gas 
model  (dashed  line)  vs.  analytical  solution  of  convective-conductive  model 
(solid  line)  in  the  limit  of  local  thermal  equilibrium.  The  cell  is  in  the  initial  stages 
of  being  heated  by  hot  air  with  velocity,  u-  1ms-1,  and  at  a  rate  K  =  1  °C  s-1 . 
Simulation  time  (in  seconds)  for  each  temperature  profile  is  indicated  as  a  label 
for  each  curve  in  the  figure. 

Several  steps  were  taken  to  validate  the  numerical  solution 
of  the  two-equation,  coupled,  solid-gas  model,  Eqs.  (8)  and 
(10).  First,  the  grid  spacing  was  reduced  until  results  were  grid 
independent.  This  requirement  was  satisfied  by  discretizing  the 
channel  length  into  100  nodes.  Second,  the  time  step  was  reduced 
until  results  were  time- step  independent,  which  occurred  for  the 
same  time-step  size  as  discussed  in  Section  2.4  for  the  CFD 
model.  Finally,  the  numerical  solution  was  validated  against  the 
analytical  solution  of  the  second  reduced-order  model  in  Section 
2.3  in  the  limit  of  local  thermal  equilibrium.  In  order  to  force 
thermal  equilibrium  between  the  gas  and  solid,  the  heat  trans¬ 
fer  coefficient  was  made  1000  times  larger  than  the  baseline, 
thermally-fully  developed  value  (x*  ->  oo)  for  laminar  flow  in  a 
duct  as  predicted  by  Eq.  (29).  This  caused  the  temperature  dif¬ 
ference  between  the  gas  and  solid  to  vanish,  and  the  resulting 
solid-gas  temperature  profile  matched  that  obtained  analytically 
by  solving  Eq.  (12).  Fig.  3  shows  results  of  this  validation  test 
for  a  cell  in  the  initial  stages  of  being  heated  by  air  flowing  at 
M=lms-1  with  the  inlet  air  temperature  increasing  at  a  rate  of 
K=  1  °C  s_  1 .  The  analytical  and  numerical  solutions  match  very 
well  for  this  case  as  well  as  several  other  cases  that  were  run, 
thus  establishing  the  validity  of  the  numerical  procedure  used 
for  integration  of  the  1st  reduced-order  model,  Eqs.  (8)  and  (10). 

Having  validated  the  numerical  scheme  for  solving  the  two- 
equation  coupled,  solid-gas  model,  the  predictions  of  heating 
time  and  maximum  temperature  gradient  were  compared  to 
the  results  of  the  CFD  Fluent  simulations.  Figs.  4  and  5  show 
heating  time,  spatial  temperature  gradient,  and  temporal  temper¬ 
ature  gradient  vs.  rate  of  inlet  temperature  rise,  K ,  for  a  heating 
air  velocity  of  u-  10ms_1.  The  predictions  of  heating  time 
and  temporal  temperature  gradient  show  excellent  agreement 


Fig.  4.  Predictions  of  heating  time  vs.  K  for  hot  air  velocity,  u-  10ms_1.  The 
solid  line  is  the  two-equation,  coupled,  solid-gas  model  predictions,  and  the 
triangles  are  data  points  obtained  from  3-D  CFD  Fluent  simulations. 

for  a  wide  range  of  K.  However,  the  model  significantly  over 
predicts  spatial  temperature  gradients  (>100%  error)  as  seen 
in  Fig.  5.  In  fact,  this  was  seen  to  be  the  case  for  most  val¬ 
ues  of  u  that  were  used  in  simulations.  Further  investigation 
revealed  that  the  excessive  temperature  gradients  occurred  near 
the  inlet  and  were  highly  sensitive  to  the  assumed  behavior  of 


Fig.  5.  Predictions  of  maximum  temperature  gradient  vs.  K  for  hot  air  veloc¬ 
ity,  u-  10ms-1.  The  solid  line  is  the  two-equation,  coupled,  solid-gas  model 
predictions,  and  the  triangles  are  results  of  3-D  CFD  Fluent  simulations.  The 
upper  plot  is  for  spatial  temperature  gradient  and  the  lower  plot  is  for  temporal 
temperature  gradient. 
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heat  transfer  coefficient,  h ,  which  increases  dramatically  near 
the  inlet.  Indeed,  the  constant- wall  temperature  correlation,  Eq. 
(29),  appeared  to  predict  significantly  larger  values  of  h  (near 
the  inlet)  than  what  was  found  by  rigorous  solution  of  the  conju¬ 
gate  mass,  momentum,  and  energy  conservation  equations  in  the 
Fluent  CFD  model.  In  reality,  h  varies  with  air  velocity,  distance 
from  the  inlet,  thermal  condition  of  the  channel  walls,  and  time, 
a  relationship  not  explicitly  known,  and  requiring  solution  of  a 
much  more  complicated  conjugate  heat  transfer  problem  (see 
for  example,  [10,1 1]).  Adding  this  significant  additional  level  of 
complexity  to  the  two-equation  coupled  solid-gas  model  would 
cut  into  its  expected  advantage  of  being  much  simpler  and  more 
computationally  efficient  than  the  3-D  Fluent  model. 

Having  established  that  the  two-equation  coupled  solid-gas 
numerical  model  is  severely  limited  in  its  ability  to  be  used  as  a 
reliable  transient  process  design  tool  for  SOFCs,  the  remainder 
of  this  paper  will  focus  on  the  analytical  model  results  and  their 
application  to  optimal  design  of  SOFC  transient  heating  process. 
These  two  analytical  reduced-order  models  -  for  purely  convec¬ 
tive  and  convective-conductive  heating  -  were  solved  for  a  range 
of  realistic  values  of  effective  rate  of  inlet  temperature  rise,  Ke ff , 
and  effective  Peclet  number,  Pe.  All  calculations  were  performed 
for  the  SOFC  geometry  shown  in  Fig.  1  and  the  material  prop¬ 
erty  values  given  in  Table  1 .  The  studied  values  of  Pe  correspond 
to  physical  velocity  in  the  range  1  <  u  <  20  m  s_1,  and  the  val¬ 
ues  of  Kq ff  covers  the  range  0.01  <  K  <  5  °C  s_1  in  terms  of  the 
actual  rate  K  of  the  inlet  air  temperature  rise.  Model  predictions 
of  dimensionless  heating  time,  and  maximum  temperature  gra¬ 
dients  (spatial  and  temporal)  are  shown  in  Figs.  6-8.  Strictly 
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Fig.  6.  Dimensionless  heating  time,  ,  as  a  function  of  Ke ff  for  various  effective 
Peclet  numbers.  The  solid  lines  are  the  reduced-order  analytical  model  predic¬ 
tions  and  the  data  points  are  results  of  3-D  Fluent  simulations.  For  Kes  1  the 
analytical  models  (regardless  of  Peclet  number)  accurately  predict  heating  time. 
As  Ke ff  increases,  the  analytical  models  are  less  accurate  but  still  predict  the 
correct  trend.  For  large  Peclet  numbers  (e.g.  Pe=  16.7)  both  analytical  models 
(i.e.,  for  purely  convective  and  convective-conductive  heating)  predict  identical 
results. 


Fig.  7.  Maximum  dimensionless  spatial  temperature  gradient  as  a  func¬ 
tion  of  Ke ff  for  various  Peclet  numbers.  The  solid  lines  are  the  analytical 
convective-conductive  model  predictions  and  the  data  points  are  results  of  3- 
D  Fluent  simulations.  Temperature  gradients  decrease  with  a  decrease  in  Pe  at 
a  fixed  Keff.  The  case  of  Pe=  16.7  is  as  calculated  from  the  purely  convective 
model  (the  limiting  case  of  Pe  — ►  oo),  which  yields  results  identical  to  those 
obtained  using  the  convective-conductive  model  at  such  a  high  Peclet  number. 

speaking,  the  purely  convective  model  is  the  limiting  case  of  the 
convective-conductive  model,  in  the  limit  of  Pe  ->  oo.  However, 
the  convective-conductive  model  results  rapidly  converge  to  the 
purely  convective  model  results  for  effective  Peclet  numbers  as 
small  as  10.  Thus,  the  results  for  the  case  of  Pe=  16.7  shown 
in  the  plots  are  obtained  using  the  purely  convective  model  cal¬ 
culations  with  the  realization  that  further  increases  in  Pe  do  not 
change  the  behavior. 


f-T0)iK 


Fig.  8.  Maximum  dimensionless  temporal  temperature  gradient  as  a  function 
of  Kq  ff  for  various  Peclet  numbers.  The  solid  lines  are  the  analytical  model 
predictions  and  the  data  points  are  results  of  3-D  Fluent  simulations.  Excellent 
agreement  between  the  analytical  models  and  3-D  results  is  observed  across  a 
wide  range  of  Keff.  At  values  of  Kcs  <  ~0.5  the  Peclet  dependence  is  negligible 
and  the  purely  convective  heating  analytical  model  provides  satisfactory  results. 
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3.1.  Heating  time 

The  first  quantity  of  interest  in  transient  thermal  modeling  of 
SOFCs  is  the  total  time  required  to  heat  the  cell  to  the  desired 
operating  temperature.  Fig.  6  shows  total  dimensionless  heat¬ 
ing  time  for  various  Keff.  The  dimensionless  heating  time  pre¬ 
dicted  by  the  purely  convective  reduced-order  model  is  simply 

=  1  / Kq ff  +  1,  which  is  a  much  ‘nicer’  result  than  the  infinite 
series  solution  of  the  convective-conductive  model.  However, 
it  is  apparent  from  Fig.  6  that  the  convective-conductive  model 
and  the  purely  convective  model  give  virtually  indistinguish¬ 
able  results  for  any  (even  small)  Pe  as  long  as  Ke ff  is  sufficiently 
small.  At  larger  KQ ff,  some  Pe  dependence  is  seen,  with  heating 
time  asymptotically  approaching  a  finite  value  (dimensionless 
value  greater  than  1)  depending  on  Pe.  Another  feature  of  Fig.  6 
is  that  excellent  agreement  (<2%  error)  is  obtained  between  the 
analytical  models  and  3-D  Fluent  simulation  results  at  small  val¬ 
ues  of  effective  rate  of  inlet  temperature  rise,  Ke ff  ^<0.1.  This 
lack  of  Pe  dependence  and  good  agreement  with  rigorous  CFD 
simulations  is  hardly  surprising  when  considered  in  the  context 
of  the  physical  time  scales,  rc  and  identified  in  Section  2.2  of 
this  paper.  Small  Ke ff  implies  that  the  time  scale  associated  with 
the  prescribed  rate  of  inlet  temperature  rise,  x\  =  (7f  —  Tq)/K,  is 
primarily  responsible  for  determining  how  fast  the  cell  can  be 
heated.  Thus,  with  respect  to  predictions  of  total  heating  time,  the 
details  of  heat  transfer  within  the  cell  are  practically  irrelevant! 
(The  relevant  portion  of  the  analysis  was  implicitly  included  in 
the  statement  that  KQ ff  is  small.)  For  large  KQ ff  this  is  no  longer 
the  case;  the  time  scale  (and  mechanism)  for  heat  flow  within 
the  channel  (whether  by  advection  or  conduction)  is  of  criti¬ 
cal  importance,  and  a  rigorous  multi-dimensional  thermal-fluid 
analysis  is  required  in  order  to  accurately  predict  heating  time. 
Fortunately,  KQ ff  (as  related  to  SOFC  unit  cells)  is  often  small, 
and  use  of  the  purely  convective  model  for  prediction  of  heating 
time  is  quite  accurate  for  practical  purposes. 

3.2.  Maximum  spatial  temperature  gradient 

The  second  quantity  of  interest  in  transient  thermal  mod¬ 
eling  of  SOFCs  is  maximum  spatial  temperature  gradient, 
which  sets  an  upper  limit  on  thermomechanical  stress  devel¬ 
oped  in  the  cell  during  its  transient  heating.  The  purely  convec¬ 
tive  model  prediction  of  dimensionless  maximum  temperature 
gradient  is  explicitly  stated  in  Eq.  (28).  On  the  other  hand, 
convective-conductive  model  predictions  of  temperature  gra¬ 
dient  were  calculated  numerically  from  the  analytical  solution 
of  Eq.  (17).  Results  of  these  two  models  along  with  predictions 
obtained  using  the  3-D  Fluent  simulations  are  shown  in  Fig.  7. 
Satisfactory  agreement  between  the  models  (<20%  error)  is  seen 
at  high  effective  Peclet  numbers  {Pe  >  8)  for  KQ ff  ~  <  1 .  It  is  clear 
that  the  purely  convective  model  over  predicts  temperature  gra¬ 
dients  for  very  low  Pe ,  although  it  still  predicts  the  correct  trends. 
The  convective-conductive  model  provides  some  improvement 
at  very  low  Pe ,  but  it  is  not  reliable  for  intermediate  values.  An 
expected,  yet  interesting  feature  of  the  graphs  shown  in  Fig.  7  is 
that  the  dimensionless  temperature  gradient  for  any  given  Ke ff 
decreases  with  a  decrease  in  Pe  from  its  maximum  at  large  Pe 


values  (purely  convective  model).  This  is  due  to  the  fact  that 
the  addition  of  conduction  as  a  significant  heat  transfer  mech¬ 
anism  reduces  temperature  gradients  along  the  axial  direction. 
This  does  not  imply,  however,  that  real  (physical)  temperature 
gradients  decrease  with  a  decrease  in  velocity.  The  only  way 
that  lowering  the  effective  Peclet  number  will  result  in  lower 
(real)  temperature  gradients  is  by  increasing  the  effective  ther¬ 
mal  diffusivity  of  the  cell.  Unfortunately,  this  may  not  be  a  viable 
design  option,  as  it  requires  changing  the  geometry  of  the  cell 
or  selecting  new  materials,  which  are  usually  optimized  for  the 
most  efficient  steady-state  operation. 

3.3.  Maximum  temporal  temperature  gradient 

One  last  quantity  of  interest  in  transient  operation  of  a  cell 
is  the  maximum  temporal  temperature  gradient  in  cell  materi¬ 
als.  This  is  important  because  the  cell  is  composed  of  layers 
made  of  different  materials  with  different  thermal  expansion 
coefficients  and  characteristic  time  scales  for  creep.  If  the  tem¬ 
poral  temperature  gradient  exceeds  the  rate  at  which  creep  can 
relax  interfacial  stress  due  to  thermal  expansion  mismatch,  then 
delamination,  cracking,  or  other  failure  may  occur.  Fig.  8  is  a 
plot  of  the  maximum  dimensionless  temporal  temperature  gra¬ 
dient  for  various  Pe  and  KQ ff .  Clearly,  for  the  values  of  KQ ff  less 
than  ~1,  the  temporal  gradient  is  essentially  independent  of  Pe , 
showing  excellent  agreement  with  the  3-D  Fluent  model  results. 
Thus,  results  given  by  the  purely  convective  heating  model  are 
valid  for  prediction  of  the  local  rate  of  temperature  change  dur¬ 
ing  transient  heating  of  the  SOFC. 

For  the  three  global  quantities  of  interest-heating  time,  maxi¬ 
mum  spatial  temperature  gradient,  and  maximum  temporal  tem¬ 
perature  gradient,  we  find  then,  that  the  purely  convective  model 
provides  satisfactory  results  within  certain  limits  of  parameters. 
For  heating  time  predictions,  the  model  is  valid  for  KQ ff  <0.1 
and  1  <Pe<o o.  The  purely  convective  model  may  be  valid 
at  even  smaller  effective  Peclet  numbers,  but  we  do  not  have 
data  to  support  this,  and  do  not  envision  scenarios  in  which 
such  low  Pe  is  a  realistic  operating  regime.  For  maximum  spa¬ 
tial  temperature  gradients,  the  model  is  reasonably  accurate  for 
Kq  ff  <  1  and  5  <  Pe  <  oo.  If  in-depth,  precise  local  information, 
beyond  simply  the  maximum  temperature  gradient,  is  required, 
we  recommend  multi-dimensional,  CFD  modeling  such  as  that 
described  in  Section  2.4.  The  purely  convective  heating  model 
is  also  valid  for  predicting  maximum  temporal  temperature 
gradients  in  the  ranges  Ke ff  <  land  l  <Pe<oo.  Because  the 
convective-conductive  model  provides  little  improvement  for 
expanding  this  range  of  validity,  we  recommend  a  design  pro¬ 
cedure  for  the  heating/cooling  process  based  on  the  purely  con¬ 
vective  model  only.  Note  that  this  model  is  shown  to  have  good 
predictive  capabilities  for  the  global  heating/cooling  character¬ 
istics  just  mentioned,  and  not  necessarily  for  giving  detailed 
time- varying  temperature  fields  within  the  cell. 

4.  Design  maps  for  transient  SOFC  heating  and  cooling 

Figs.  6-8  are  design  maps,  which  can  be  used  to  develop 
a  protocol  for  an  optimal  SOFC  heating  process.  Typically, 
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the  cell  geometry  and  material  selection  will  be  optimized  for 
steady-state  cell  performance,  leaving  the  hot  air  velocity,  u ,  and 
heating  rate,  K ,  as  control  variables  for  optimizing  the  transient 
process.  Assuming  that  the  maximum  allowable  stress  (and  cor¬ 
responding  temperature  gradient)  has  been  specified  through  a 
thermo-mechanical  stress  analysis  of  the  unit  cell,  the  design 
goal  is  to  minimize  heating  time  through  proper  choice  of  K 
and  u.  Because  of  the  use  of  properly  scaled  dimensionless  vari¬ 
ables,  Figs.  6-8  have  the  advantage  of  presenting  a  large  amount 
of  information  on  a  single  curve.  A  detailed,  step-by-step  design 
procedure  based  on  similar  maps  has  been  presented  by  Damm 
and  Fedorov  [7] .  Here,  however,  it  is  advantageous  to  refer  back 
to  physical  parameters  in  outlining  a  design  procedure  for  opti¬ 
mization  of  the  heating/cooling  process. 

In  particular,  Fig.  9  shows  a  design  map  created  from  the 
purely  convective  model  predictions  (see  Eqs.  (27)  and  (28)). 
This  model  was  shown  in  Section  3  to  be  valid  over  specific 
ranges  of  effective  Peclet  number  and  effective  rate  of  inlet  tem¬ 
perature  rise.  The  map  in  Fig.  9  is  specific  to  the  geometry  and 
material  selection  previously  specified  (see  Fig.  1  and  Table  1), 
and  considers  a  unit  cell  being  heated  from  25  to  625  °C  (the 
map  is  applicable  for  any  process  involving  a  600  °C  tempera¬ 
ture  change).  Lines  of  constant  K  and  u  are  shown  with  maximum 
spatial  temperature  gradient  on  the  vertical  axis  and  heating  time 
on  the  horizontal  axis.  Once  a  horizontal  line  corresponding  to 
the  maximum  allowable  temperature  gradient  is  specified,  K  and 
u  can  be  chosen  below  that  line  such  that  heating  time  is  mini¬ 
mized.  This  allows  some  flexibility  in  choice  of  K  and  u  but  these 
parameters  are  likely  to  be  restricted  by  system  constraints  such 
as  pumping  power,  heater  size,  etc.  Alternatively,  the  maximum 
allowable  heating  time  could  be  the  specified  design  constraint, 
in  which  case  K  and  u  are  chosen  such  that  temperature  gradient 
is  minimized.  This  would  give  an  estimate  of  the  allowable  gra¬ 
dients  that  the  cell  must  be  designed  to  withstand  if  the  heating 


Fig.  9.  Design  map  based  on  purely  convective  model  predictions  of  spatial 
temperature  gradient  (°Cm-1)  vs.  heating  time  (s)  for  various  K  and  u.  This 
map  is  specific  to  the  cell  geometry  and  materials  described  in  Fig.  1  and  Table  1, 
and  is  applicable  to  a  cell  being  heated  or  cooled  over  a  600  °C  temperature 
interval. 


time  requirement  is  to  be  met.  Notice  that  increasing  the  air  flow 
speed,  u ,  always  allows  movement  in  a  favorable  direction  on  the 
map  (i.e.,  reduction  in  both  the  heating  time  and  maximum  tem¬ 
perature  gradients),  while  increasing  K  yields  the  mixed  result 
of  lowered  heating  time  but  increased  temperature  gradients.  A 
design  map  such  as  Fig.  9  is  thus  an  excellent  graphical  repre¬ 
sentation  of  how  the  design  space  may  be  efficiently  searched 
to  yield  approximate  values  of  hot  air  velocity  and  heating  rate 
to  achieve  an  optimal  design  of  SOFC  transient  heating  process. 
Fig.  9  is  equally  valid  for  heating  and  cooling  of  the  cell;  how¬ 
ever,  the  maximum  allowable  temperature  gradient  for  a  heating 
process  may  be  different  than  that  for  cell  cooling  due  to  differ¬ 
ences  in  thermomechanical  material  behavior  under  heating  and 
cooling  conditions.  Once  the  parameter  K  and  u  values  result¬ 
ing  in  optimal  transient  operation  have  been  obtained,  only  a 
minimal  amount  of  highly-intensive  CFD  computations  can  be 
used  to  yield  the  detailed  information  required  to  complete  a 
thermo-mechanical  failure  analysis  of  the  cell. 

5.  Conclusions  and  future  work 

Reduced-order  transient  thermal  models  of  varying  complex¬ 
ity  have  been  considered  for  optimizing  the  heating  and  cooling 
of  an  SOFC  unit  cell.  The  first  and  most  general  one,  the  two- 
equation,  coupled  solid-gas  model  resulted  in  fast  numerical 
solution  of  the  problem,  compared  to  CFD  simulations,  but  was 
shown  to  be  unreliable  for  predicting  maximum  temperature  gra¬ 
dients.  The  next  two  models  the  convective-conductive  heating 
and  purely  convective  heating  permitted  analytical  solutions  of 
the  time- varying  temperature  field  in  the  cell.  The  simplest  of 
these,  the  purely  convective  heating  model,  also  yielded  explicit, 
algebraic  relationships  between  heating  time,  temperature  gra¬ 
dients,  hot  air  velocity,  and  heating  rate.  These  predictions  of 
integral  thermal  quantities  (i.e.,  the  heating  time  and  maximum 
spatial  and  temporal  temperature  gradients)  were  shown  to  be 
valid  for  a  fairly  broad  range  of  operating  parameters  through 
comparison  with  fully  3-D  Fluent  simulations.  The  more  general 
convective-conductive  model  of  the  cell  provided  little  improve¬ 
ment  in  accuracy  over  the  purely  convective  model;  however,  its 
formulation  and  analysis  led  to  the  definition  of  the  appropriate 
scales  for  physical  velocity  and  the  dimensionless  parameters 
that  govern  cell  transient  response. 

Predictions  of  the  analytical  models  have  been  presented  in 
the  form  of  generalized  thermal  design  maps,  and  a  specific 
example  of  a  design  map  based  on  the  purely  convective  model 
was  used  to  develop  a  conceptual  framework  for  optimizing  a 
heating  and  cooling  process.  Our  analysis  shows  that  increasing 
the  velocity  of  the  hot  air  stream,  and  lowering  the  Peclet  number 
(by  increasing  the  effective  thermal  diffusivity  of  the  cell)  leads 
to  the  optimal  design,  which  minimizes  heating  time  under  the 
constraint  of  maximum  allowable  temperature  gradients.  While 
this  result  is  hardly  surprising  in  hindsight,  the  ability  of  the 
purely  convective  model  to  accurately  predict  favorable  design 
trends  with  little  or  no  computational  effort  makes  it  a  powerful 
tool  for  searching  the  design  space  in  the  early  stages  of  transient 
process  development. 
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Throughout  the  development  of  these  reduced-order  models 
and  thermal  design  maps,  it  has  been  assumed  that  the  maximum 
allowable  temperature  gradient  (to  avoid  failure  of  cell  com¬ 
ponents)  is  known  a  priori  based  on  thermo-mechanical  stress 
analysis.  This  is  not  necessarily  the  case  and  efforts  similar  to 
[1]  are  required  for  the  development  of  relationships  between 
failure,  maximum  allowable  stresses,  and  maximum  allowable 
temperature  gradients  as  relevant  to  SOFC  transients.  Also,  the 
effects  of  thermal  cycling  require  additional  research,  in  particu¬ 
lar  quantification  of  cell  lifetime  in  terms  of  number  of  cycles  and 
maximum  temperature  gradients  developed  during  each  cycle. 

Another  assumption  of  the  current  analysis  that  will  be  chal¬ 
lenged  is  that  the  cell  is  perfectly  insulated.  In  reality,  the  cell 
heat  losses  are  temperature  dependent  (increasing  as  the  cell 
heats  up),  which  raises  the  possibility  that  the  heating  times 
could  be  much  longer  than  what  has  been  predicted  here.  In  the 
most  extreme  case,  a  cell  might  never  reach  the  desired  operat¬ 
ing  temperature  if  the  magnitude  of  heat  input  is  not  sufficiently 
larger  than  heat  losses.  The  assumption  of  negligible  flow  in 
the  fuel  channel  may  also  come  into  question  depending  on  the 
flow  rate  of  a  reducing  gas  to  keep  the  anode  from  oxidation 
at  elevated  temperatures  (i.e.  T>  400  °C).  To  account  for  this, 
the  energy  conservation  equation  for  the  fuel  channel  in  the  sys¬ 
tem  of  Eqs.  (5)  and  associated  boundary/inlet  condition  would 
have  to  be  modified  to  include  an  advective  term.  As  a  result, 
the  two-equation  model  (Section  2.1)  would  become  a  three- 
equation  model.  This  extension  is  straight-forward  and  follows 
the  general  theoretical  framework  we  have  outlined.  Alterna¬ 
tively,  the  heat  losses  due  to  the  flowing  reducing  gas  on  the 
anode  side  can  be  included  into  the  two-equation  and  single 
equation  models  by  using  an  ad-hoc  approach  through  specify¬ 
ing  a  semi-empirical  heat  loss  term  along  the  unit  cell  length.  In 
either  case,  the  treatment  is  highly  case-specific,  so  we  refrain 
from  its  detailed  consideration  in  presenting  the  results  of  a 
more  general  nature.  Finally,  extension  and  validation  of  the  cur¬ 
rent  analysis  for  other  geometries,  such  as  cross-flow  cells,  will 
greatly  expand  the  applicability  of  the  results.  This  future  work 
will  complement  and  enhance  the  presently  developed  simple 
and  computationally  efficient  design  rules  for  transient  heat¬ 
ing/cooling  of  the  SOFC. 
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Look  for  a  solution  of  the  form 
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where  A  and  B  are  arbitrary  constants.  It  can  be  shown  that  if 
A  =  Pe/2  and  B  =  Pel  A  then  the  advection  term  in  Eq.  (30)  van¬ 
ishes,  yielding 
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where  the  dimensionless  inlet  forcing  function  from  Eq.  (31) 
becomes 
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The  new  problem  can  be  expressed  as 
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where  the  unknown  functions  are  selected  such  that  the  boundary 
conditions  become  homogeneous  yielding 
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which  can  now  be  solved  analytically.  The  homogeneous  por¬ 
tion  of  Eq.  (37)  is  solved  using  separation  of  variables  and  the 
complete  solution  is  assumed  a  Fourier  series  expansion  based 
on  the  complete  eigenvector  set  of  the  homogeneous  problem 
but  with  time  varying  coefficients  Cn(t) 


Appendix  A.  Analytical  solution  of 
convective-conductive  Eq.  (17) 
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We  seek  a  solution  of  the  dimensionless  governing  Eq.  (17) 
(*  superscript  has  been  dropped  for  clarity) 
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The  denominator  is  the  squared  eigenvector  norm,  which  we 
will  call  Gn  for  simplicity 


Finally,  Cn(t)  for  the  final,  constant  portion  of  the  inlet  temper¬ 
ature  function  (see  Eq.  (1 1))  is  found  to  be 
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where  t\  =  l/^eff  is  the  simulation  time  after  which  the  inlet 
temperature  is  held  constant. 


with  the  eigenvalues  con  =  Xn  \J~Pe  satisfying  the  following  char¬ 
acteristic  equation: 
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It  follows,  after  some  manipulation  [7]: 
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After  applying  the  inlet  temperature  function  from  Eq.  (1 1)  and 
carrying  out  the  integration,  Cn(t)  is 
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